TetMesh Geodesics

import numpy as np
from lapy import TetMesh
from lapy.plot import plot_tet_mesh
import plotly
# plotly.offline.init_notebook_mode(connected=True)
import plotly.io as pio
pio.renderers.default = "sphinx_gallery"

T = TetMesh.read_vtk('../data/cubeTetra.vtk')
#T.is_oriented()
T.orient_()

--> VTK format         ...
 --> DONE ( V: 9261 , T: 48000 )

Flipped 24000 tetrahedra

24000

from lapy import Solver
fem = Solver(T,lump=True)

evals, evec = fem.eigs(10)

TetMesh with regular Laplace
Solver: spsolve (LU decomposition) ...

# also get A,B (lumped), and inverse of B (easy as it is diagonal)
A, B = fem.stiffness, fem.mass
Bi = B.copy()
Bi.data **= -1

evnum=1
cutting = ['x<0.5']
# also here we comment all plots to reduce file size
# uncomment and take a look
plot_tet_mesh(T,vfunc=evals[evnum]*evec[:,evnum],plot_edges=None,plot_levels=False,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)

Mesh is oriented, nothing to do
Found 3040 triangles on boundary.

from lapy.diffgeo import compute_gradient
from lapy.diffgeo import compute_divergence
grad = compute_gradient(T,evec[:,evnum])
divx = -compute_divergence(T,grad)
vfunc = Bi*divx

cutting = ['x<0.5']
plot_tet_mesh(T,vfunc=vfunc,plot_edges=None,plot_levels=False,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)

Mesh is oriented, nothing to do
Found 3040 triangles on boundary.

np.max(np.abs(vfunc-(evals[evnum]*evec[:,evnum])))
dd = np.abs(vfunc-(evals[evnum]*evec[:,evnum]))
max(dd)

0.006755829

from lapy import heat

tria = T.boundary_tria()
bvert = np.unique(tria.t)

u = heat.diffusion(T,bvert,m=1)
cutting = ['x<0.5']
plot_tet_mesh(T,vfunc=u,plot_edges=None,plot_levels=True,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)

Found 4800 triangles on boundary.
TetMesh with regular Laplace
Matrix Format now:  csc
Solver: spsolve (LU decomposition) ...
Mesh is oriented, nothing to do
Found 3040 triangles on boundary.

# compute gradient of heat diffusion, normalize it, and compute the divergence of normalized gradient
tfunc = compute_gradient(T, u)

# flip and normalize it
X = -tfunc / np.sqrt((tfunc ** 2).sum(1))[:, np.newaxis]
X = np.nan_to_num(X)

# compute divergence
divx = compute_divergence(T, X)

# compute distance
from scipy.sparse.linalg import splu
useCholmod = True
try:
    from sksparse.cholmod import cholesky
except ImportError:
    useCholmod = False

A, B = fem.stiffness, fem.mass # computed above when creating Solver

H=A
b0=-divx

# solve H x = b0
print("Matrix Format now: "+H.getformat())
if useCholmod:
    print("Solver: cholesky decomp - performance optimal ...")
    chol = cholesky(H)
    x = chol(b0)
else:
    print("Solver: spsolve (LU decomp) - performance not optimal ...")
    lu = splu(H)
    x = lu.solve(b0)

x = x - np.min(x)

Matrix Format now: csc
Solver: spsolve (LU decomp) - performance not optimal ...

cutting = ['x<0.5']
plot_tet_mesh(T,vfunc=x,plot_edges=None,plot_levels=True,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)
max(x), 0.5*np.sqrt(3.0)

Mesh is oriented, nothing to do
Found 3040 triangles on boundary.

(0.69931483, 0.8660254037844386)

#debug gradient
v1func =  T.v[:,0]* T.v[:,0] + T.v[:,1]* T.v[:,1] + T.v[:,2]* T.v[:,2]# x coord

grad = compute_gradient(T,v1func)
glength = np.sqrt(np.sum(grad * grad, axis=1))
fcols=glength
fcols=grad[:,2]
cutting = ['x<0.5']
cutting = None
plot_tet_mesh(T,vfunc=None,tfunc=fcols,plot_edges=None,plot_levels=False,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False)

Mesh is oriented, nothing to do
Found 4800 triangles on boundary.

divx = compute_divergence(T, grad)
divx2 = Bi * divx
cutting = ['z<0.5']
plot_tet_mesh(T,vfunc=divx2,plot_edges=True,plot_levels=False,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)

Mesh is oriented, nothing to do
Found 3040 triangles on boundary.

divx2[5000:5010]

array([5.9999948, 6.0000215, 6.0000215, 5.999988 , 6.000053 , 5.999975 ,
       5.9999676, 6.000024 , 6.000013 , 6.000008 ], dtype=float32)

T.avg_edge_length()

0.06365621

from lapy.diffgeo import compute_geodesic_f
from lapy import heat

tria = T.boundary_tria()
bvert=np.unique(tria.t)

# get heat diffusion
u = heat.diffusion(T,bvert, m=1)

gu = compute_geodesic_f(T,u)

cutting = ['x<0.5']
plot_tet_mesh(T,vfunc=gu,plot_edges=None,plot_levels=True,cutting=cutting,edge_color='rgb(50,50,50)',html_output=False,flatshading=True)

Found 4800 triangles on boundary.
TetMesh with regular Laplace
Matrix Format now:  csc
Solver: spsolve (LU decomposition) ...
TetMesh with regular Laplace
Matrix Format now: csc
Solver: spsolve (LU decomposition) ...
Mesh is oriented, nothing to do
Found 3040 triangles on boundary.